Homework 3

Name:

Collaborators:

Homework due Friday 10/23/2015 before 23:55

Please, write the name of your collaborators. We encourage collaboration for solving the homework; however, you need to write your own code and acknowledge your collaborators.

Question 1: Mueller's Method

Mueller's method is another root-finding method. It uses 3 points to approximate your function locally with a parabola, and then computes the intersection of the parabola with $y=0$ to update one of your points. The main advantage of Mueller's method with respect to Newton's method and its variant, is that it can compute complex roots.

1.a You will write a function rootMuellerMethod with inputs: $f$ a function handle, $p0$, $p1$, $p0$ the initial guesses, $\epsilon$ the tolerance, and $Nmax$ the maximum number of iterations.

Your function will raise an error if the convergence is not reached in the specified number of iterations.

The output of the function will depend on the optional parameter $history$. If $history$ is false, the function will output the last guess $p$; on the other hand, if $history$ is true the output will be vector $p$ with all the intermediate steps.

Hint: you may need to use complex arithmetic. In order to avoid errors, you need to force p0,p1 and p2 to be complex numbers.


In [1]:
function rootMuellerMethod(f,p0,p1,p2, ϵ, Nmax; history = false )
    p0, p1, p2 = (p0 + 0*im,p1 + 0*im,p2 + 0*im)  # forcing p0,p1,p2 to be complex numbers
    history && (hist = ([p0 p1 p2][:]))  # defining the history vector
    # and adding the first 3 guesses
    # the code is extract straight from the book
    h1 = p1 - p0
    h2 = p2 - p1
    δ1 = (f(p1) - f(p0))/h1
    δ2 = (f(p2) - f(p1))/h2
    d  = (δ2-δ1)/(h2 + h1)
    
    for i = 1:Nmax
        b = δ2 + h2*d
        # compute the deternminant
        # we need to use complex arithmetic here
        D = sqrt(b^2 - 4*f(p2)*d)
        # we select the appropiate root of the parabola
        abs(b-D) < abs(d+D) ? E = b+D : E = b-D
        h = -2*f(p2)/E
        p = p2+h
        history && (push!(hist,p))
        # if the update is smaller than the tolerance
        # we give back the answer
        if (abs(h) < ϵ)
            # if history is true we return the vector
            # with all the intermediate steps; otherwise
            # we retunr just the last p 
            history? (return hist) : (return p)
        end
        # we get ready for the next iteration
        p0 = p1
        p1 = p2
        p2 = p
        h1 = p1 - p0
        h2 = p2 - p1
        δ1 = (f(p1) - f(p0))/h1
        δ2 = (f(p2) - f(p1))/h2
        d  = (δ2-δ1)/(h2 + h1)
    end
    error("Method failed after Max Number of iterations")
end


Out[1]:
rootMuellerMethod (generic function with 1 method)

Question 2: Comparison of the Root-Finding algorithms

2.a You will write a script to find all the roots to within $10^{-6}$ of the following polynomials using Mueller's Method. You will use the function that you just wrote, with suitable initial guesses. If you can not find all the complex roots you may want to take a look at Theorem 2.20 in your textbook. In Julia you can use the conj function to compute the conjugate of a complex number.


In [2]:
p1(x) = x.^3 - 2*x.^2 - 5
p2(x) = x.^3 + 3*x.^2 - 1
p3(x) = x.^3 - x - 1
p4(x) = x.^4 + x.^2 - x - 3
p5(x) = x.^4 + 4.001*x.^2 + 4.002*x + 1.101
p6(x) = x.^5 - x.^4 + 2*x.^2 + x - 4


Out[2]:
p6 (generic function with 1 method)

In order to make the correction of the Homework easier, you will print in the console (using the function print) the roots for each polynomial.


In [11]:
# write your script here
# roots of p1
root1_1 = rootMuellerMethod(p1,0.1,-1.0, 3.0, 1e-6, 30, history = true )
root1_2 = rootMuellerMethod(p1,0.1,-1.0,-3.0, 1e-6, 30, history = true )
print("Roots of P1 \n ",root1_1[end],  "\n")
print(root1_2[end],  "\n")
print(root1_2[end]',  "\n")
# roots of p2
root1_1 = rootMuellerMethod(p2,0.1,-1.0, 3.0, 1e-6, 30, history = true )
root1_2 = rootMuellerMethod(p2,0.1,-1.0,-3.0, 1e-6, 30, history = true )
root1_3 = rootMuellerMethod(p2,1,-1.0,-3.0, 1e-6, 30, history = true )
print("Roots of P2 \n ",root1_1[end],  "\n")
print(root1_2[end],  "\n")
print(root1_3[end],  "\n")
# roots of p3
root1_1 = rootMuellerMethod(p3,0.1,-1.0,-3.0, 1e-6, 30, history = true )
root1_2 = rootMuellerMethod(p3,1,-1.0,-3.0, 1e-6, 30, history = true )
print("Roots of P3 \n",root1_1[end],  "\n")
print(root1_1[end]',  "\n")
print(root1_2[end],  "\n")

# roots of p4
root1_1 = rootMuellerMethod(p4,-11,1.0, 20.0, 1e-6, 30, history = true )
root1_2 = rootMuellerMethod(p4,0.1,-1.0,-3.0, 1e-6, 30, history = true )
root1_3 = rootMuellerMethod(p4,-2.3,-2.5,-1.0, 1e-6, 200, history = true )
print("Roots of P4 \n",root1_1[end],  "\n")
print(root1_1[end]',  "\n")
print(root1_2[end],  "\n")
print(root1_3[end],  "\n")

# roots of p5
root1_1 = rootMuellerMethod(p5,-11,1.0, 20.0, 1e-6, 30, history = true )
root1_2 = rootMuellerMethod(p5,0.1,-1.0,-3.0, 1e-6, 30, history = true )
print("Roots of P5 \n",root1_1[end],  "\n")
print(root1_1[end]',  "\n")
print(root1_2[end],  "\n")
print(root1_2[end]',  "\n")
# roots of p6
root1_1 = rootMuellerMethod(p6,1.1,1.0, 1.2, 1e-6, 30, history = true )
root1_2 = rootMuellerMethod(p6,0.8,0.85, .9, 1e-6, 30 ,history = true )
root1_3 = rootMuellerMethod(p6,-2.3,-2.5,-1.0, 1e-6, 200, history = true )
print("Roots of P6 \n",root1_1[end],  "\n")
print(root1_2[end],  "\n")
print(root1_2[end],  "\n")
print(root1_3[end],  "\n")


Roots of P1 
 2.6906474480286136 + 0.0im
-0.3453237240142915 + 1.3187267795713722im
-0.3453237240142915 - 1.3187267795713722im
Roots of P2 
 0.5320888862379617 - 2.4635870169018266e-14im
-2.879385241571816 + 0.0im
-0.65270364466614 + 0.0im
Roots of P3 
-0.6623589786226283 - 0.5622795120620463im
-0.6623589786226283 + 0.5622795120620463im
1.3247179572439298 - 4.4693717316438534e-13im
Roots of P4 
-0.13784110182551523 - 1.527312250886601im
-0.13784110182551523 + 1.527312250886601im
1.2756822036508013 + 0.0im
-1.0 + 0.0im
Roots of P5 
0.45855586303572055 + 2.095858592591658im
0.45855586303572055 - 2.095858592591658im
-0.4585558630371269 - 0.17006974716965353im
-0.4585558630371269 + 0.17006974716965353im
Roots of P6 
1.1430034814355925 + 0.0im
1.1430034814355925 + 0.0im
1.1430034814355925 + 0.0im
1.1430034814355925 - 1.049200612163558e-16im

2.c You will write a script that finds all the roots of the polynomials, using the Newton's Method. You will use the Newton's method you wrote in Question 2 from Homework 2. Can you find all roots?


In [4]:
function newtonMethod(f,dfdx, p0, ϵ, Nmax; history=false)
    # if history true, declare an array an push p0 inside
    history && (pHist = [p0][:])
    # for loop
    for i = 1:Nmax
        p = p0 - f(p0)/dfdx(p0) # p_{n+1} = p_n - f(p_n)/f'(p_n)
        history && push!(pHist, p) # if history is true, push p_{n+1} into the pHist
        if abs(p-p0)<ϵ  # if the update is smaller than the tolerance return the answer
            history? (return pHist): (return p)
        end
        p0 = p # update p0 to start a new iteration
    end
    error("Accuracy not achieved within the specified number of iterations")
end


Out[4]:
newtonMethod (generic function with 1 method)

In order to make the correction of the Homework easier, you will print in the console (using the function print) the roots for each polynomial. You need to define by hand the derivative of the polynomial in order to be passed as a function handle to your Newton's method.


In [10]:
# write your script here
# roots of p1
dp1dx(x) = 3x.^2 - 4*x 
rootNewton1 = newtonMethod(p1,dp1dx, 2.0, 1e-6, 30)
println("Reals roots of P1 \n", rootNewton1)
# roots of p2
dp2dx(x) = 3x.^2 + 6*x 
rootNewton2 = newtonMethod(p2,dp2dx, 2.0, 1e-6, 30)
println("Reals roots of P2 \n", rootNewton2)
# roots of p3
dp3dx(x) = 3x.^2 - 1 
rootNewton3 = newtonMethod(p3,dp3dx, 2.0, 1e-6, 30)
println("Reals roots of P3 \n", rootNewton3)

# roots of p4
dp4dx(x) = 4x.^3 + 2*x -1 
rootNewton4 = newtonMethod(p4,dp4dx, 2.0, 1e-6, 30)
rootNewton4v2 = newtonMethod(p4,dp4dx, -2.0, 1e-6, 30)
println("Reals roots of P4 \n", rootNewton4)
println( rootNewton4v2)
# roots of p6
dp6dx(x) = 5x.^4 - 4*x^3 + 4*x +1 
rootNewton6 = newtonMethod(p6,dp6dx, 2.0, 1e-6, 30)
println("Reals roots of P6 \n", rootNewton6)


Reals roots of P1 
2.6906474480286136
Reals roots of P2 
0.532088886237956
Reals roots of P3 
1.324717957244746
Reals roots of P4 
1.275682203650985
-1.0000000000000606
Reals roots of P6 
1.1430034814355925

Answer: The implementation of the Newton's method doesn't utilize complex arithmetic. In particular, we can only use real initial guesses, which means that we can only find real roots.

2.c You will write a script that finds all the roots of the polynomials, using the Companion matrix method. You will use the companion matrix algorithm you wrote in Question 3 from Homework 2.


In [6]:
# use this space to write your companion matrix root-finding method
function companionMatrix(alpha)
    n = length(alpha)-1  # the deegre of the polynomial
    M = diagm(ones(n-1,1)[:],-1) # we build and fill the matrix with ones in the first subdiagonal
    M[:,end] = -alpha[1:end-1]/alpha[end]  # we put the coefficient in the righ-most column of M
    return M
end

function rootFindingCompanionMatrix(alpha)
    (root, eigVec) = eig(companionMatrix(alpha)) 
    return root
end


Out[6]:
rootFindingCompanionMatrix (generic function with 1 method)

In order to make the correction of the Homework easier, you will print in the console (using the function print) the roots for each polynomial.


In [8]:
# write your script ro find the roots in here
# roots of p1
roots1 = rootFindingCompanionMatrix([-5 0 -2 1])
println("Roots of p1 \n",roots1)
# roots of p2
roots2 = rootFindingCompanionMatrix([-1 0 3 1])
println("Roots of p2 \n",roots2)
# roots of p3
roots3 = rootFindingCompanionMatrix([-1 -1 0 1])
println("Roots of p3 \n",roots3)
# roots of p4
roots4 = rootFindingCompanionMatrix([-3  -1 1 0 1])
println("Roots of p4 \n",roots4)
# roots of p5
roots5 = rootFindingCompanionMatrix([1.101 4.002 4.001 0 1])
println("Roots of p5 \n",roots5)
# roots of p6
roots6 = rootFindingCompanionMatrix([-4 1 2  0 -1 1 ])
println("Roots of p6 \n",roots6)


Roots of p1 
Complex{Float64}[-0.34532372401430633 + 1.3187267795713242im,-0.34532372401430633 - 1.3187267795713242im,2.6906474480286136 + 0.0im]
Roots of p2 
[0.5320888862379564,-0.6527036446661387,-2.8793852415718164]
Roots of p3 
Complex{Float64}[-0.6623589786223731 + 0.5622795120623013im,-0.6623589786223731 - 0.5622795120623013im,1.324717957244746 + 0.0im]
Roots of p4 
Complex{Float64}[1.2756822036509847 + 0.0im,-0.9999999999999998 + 0.0im,-0.13784110182549236 + 1.52731225088663im,-0.13784110182549236 - 1.52731225088663im]
Roots of p5 
Complex{Float64}[-0.4585558630371269 + 0.1700697471696595im,-0.4585558630371269 - 0.1700697471696595im,0.4585558630371271 + 2.095858592593759im,0.4585558630371271 - 2.095858592593759im]
Roots of p6 
Complex{Float64}[-0.9991033032813696 + 0.6643596348511132im,-0.9991033032813696 - 0.6643596348511132im,0.9276015625635745 + 1.2531986372341644im,0.9276015625635745 - 1.2531986372341644im,1.1430034814355923 + 0.0im]

2.d Compare the algorithms (Mueller's method, Newton's Method, and Companion matrix) in this three aspects:

  • memory footprint (how much memory is being allocated) with respect to the degree of the polynomial

  • complexity (how many operations are performed) with respect to the degree of the polynomail

  • implementation, which one is easier to implement and run.

Hints:

  • The complexity of the eigenvalue decomposition is $\mathcal{O}(p^3)$ for a $p \times p$ matrix.
  • You can find useful to use "big O" notation. $\mathcal{O}(1)$: independent of n, $\mathcal{O}(n)$: linear in n, $\mathcal{O}(n^2)$: quadratic in n, etc.
  • You can make some assumptions with respect to the number of iterations, to simplify your claims.

Answer: We denote the degree of the polynomial as $n$.

Memory:

  • Newton's Method: each call needs $\mathcal{O}(1)$ memory, but we need to store $n$ roots, then the memory footprint is $\mathcal{O}(n)$;
  • Mueller's Method: same as Newton's Method;
  • Companion Matrix: we need to create a $n\times n$ matrices, which will use $\mathcal{O}(n^2)$ memory. (if we use sparse matrices the memory can go down to $\mathcal{O}(n)$

Complexity:

We suppose that each call to $f$ and $f'$ can be done in a time independent of the degree of the polynomial $n$, i.e. in $\mathcal{O}(1)$ time; moreover, we assume that the number of iterations for convergence of the Newton and Mueller's method do not depend on the degree of the polynomail. Under these assumptions we have that:

  • Newton's Method: $\mathcal{O}(1)$ per call, but given that the polynomial has $n$ roots we need $n$ calls, then the complexity is $\mathcal{O}(n)$;
  • Mueller's Method: same as Newton's Method;
  • Companion Matrix: We need to compute the eigenvalues of a $n\times n$ matrices, using a Hessenberg reduction and the QR algorithm (you will study it later in class) we have that the complexity is $\mathcal{O}(n^3)$.

Implementation:

Newton's method can not obtain complex roots starting from a real guess. We would need to change the implementation in order to use Complex arithmetic. You can obtain all the roots with Mueller's method, even the complex ones, but you need to find suitable initial guesses, which can be labor intesive. The companion matrix is the de facto method to find roots, it is extremely easy to code and you obtain all the roots without any user intervention. However, its asymptotic cost is higher.